home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zsteqr.f < prev    next >
Text File  |  1997-06-25  |  15KB  |  505 lines
  1.       SUBROUTINE ZSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
  2. *
  3. *  -- LAPACK routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       CHARACTER          COMPZ
  10.       INTEGER            INFO, LDZ, N
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
  14.       COMPLEX*16         Z( LDZ, * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  ZSTEQR computes all eigenvalues and, optionally, eigenvectors of a
  21. *  symmetric tridiagonal matrix using the implicit QL or QR method.
  22. *  The eigenvectors of a full or band complex Hermitian matrix can also
  23. *  be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this
  24. *  matrix to tridiagonal form.
  25. *
  26. *  Arguments
  27. *  =========
  28. *
  29. *  COMPZ   (input) CHARACTER*1
  30. *          = 'N':  Compute eigenvalues only.
  31. *          = 'V':  Compute eigenvalues and eigenvectors of the original
  32. *                  Hermitian matrix.  On entry, Z must contain the
  33. *                  unitary matrix used to reduce the original matrix
  34. *                  to tridiagonal form.
  35. *          = 'I':  Compute eigenvalues and eigenvectors of the
  36. *                  tridiagonal matrix.  Z is initialized to the identity
  37. *                  matrix.
  38. *
  39. *  N       (input) INTEGER
  40. *          The order of the matrix.  N >= 0.
  41. *
  42. *  D       (input/output) DOUBLE PRECISION array, dimension (N)
  43. *          On entry, the diagonal elements of the tridiagonal matrix.
  44. *          On exit, if INFO = 0, the eigenvalues in ascending order.
  45. *
  46. *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
  47. *          On entry, the (n-1) subdiagonal elements of the tridiagonal
  48. *          matrix.
  49. *          On exit, E has been destroyed.
  50. *
  51. *  Z       (input/output) COMPLEX*16 array, dimension (LDZ, N)
  52. *          On entry, if  COMPZ = 'V', then Z contains the unitary
  53. *          matrix used in the reduction to tridiagonal form.
  54. *          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
  55. *          orthonormal eigenvectors of the original Hermitian matrix,
  56. *          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
  57. *          of the symmetric tridiagonal matrix.
  58. *          If COMPZ = 'N', then Z is not referenced.
  59. *
  60. *  LDZ     (input) INTEGER
  61. *          The leading dimension of the array Z.  LDZ >= 1, and if
  62. *          eigenvectors are desired, then  LDZ >= max(1,N).
  63. *
  64. *  WORK    (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
  65. *          If COMPZ = 'N', then WORK is not referenced.
  66. *
  67. *  INFO    (output) INTEGER
  68. *          = 0:  successful exit
  69. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  70. *          > 0:  the algorithm has failed to find all the eigenvalues in
  71. *                a total of 30*N iterations; if INFO = i, then i
  72. *                elements of E have not converged to zero; on exit, D
  73. *                and E contain the elements of a symmetric tridiagonal
  74. *                matrix which is unitarily similar to the original
  75. *                matrix.
  76. *
  77. *  =====================================================================
  78. *
  79. *     .. Parameters ..
  80.       DOUBLE PRECISION   ZERO, ONE, TWO, THREE
  81.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
  82.      $                   THREE = 3.0D0 )
  83.       COMPLEX*16         CZERO, CONE
  84.       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
  85.      $                   CONE = ( 1.0D0, 0.0D0 ) )
  86.       INTEGER            MAXIT
  87.       PARAMETER          ( MAXIT = 30 )
  88. *     ..
  89. *     .. Local Scalars ..
  90.       INTEGER            I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
  91.      $                   LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
  92.      $                   NM1, NMAXIT
  93.       DOUBLE PRECISION   ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
  94.      $                   S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
  95. *     ..
  96. *     .. External Functions ..
  97.       LOGICAL            LSAME
  98.       DOUBLE PRECISION   DLAMCH, DLANST, DLAPY2
  99.       EXTERNAL           LSAME, DLAMCH, DLANST, DLAPY2
  100. *     ..
  101. *     .. External Subroutines ..
  102.       EXTERNAL           DLAE2, DLAEV2, DLARTG, DLASCL, DLASRT, XERBLA,
  103.      $                   ZLASET, ZLASR, ZSWAP
  104. *     ..
  105. *     .. Intrinsic Functions ..
  106.       INTRINSIC          ABS, MAX, SIGN, SQRT
  107. *     ..
  108. *     .. Executable Statements ..
  109. *
  110. *     Test the input parameters.
  111. *
  112.       INFO = 0
  113. *
  114.       IF( LSAME( COMPZ, 'N' ) ) THEN
  115.          ICOMPZ = 0
  116.       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
  117.          ICOMPZ = 1
  118.       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
  119.          ICOMPZ = 2
  120.       ELSE
  121.          ICOMPZ = -1
  122.       END IF
  123.       IF( ICOMPZ.LT.0 ) THEN
  124.          INFO = -1
  125.       ELSE IF( N.LT.0 ) THEN
  126.          INFO = -2
  127.       ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
  128.      $         N ) ) ) THEN
  129.          INFO = -6
  130.       END IF
  131.       IF( INFO.NE.0 ) THEN
  132.          CALL XERBLA( 'ZSTEQR', -INFO )
  133.          RETURN
  134.       END IF
  135. *
  136. *     Quick return if possible
  137. *
  138.       IF( N.EQ.0 )
  139.      $   RETURN
  140. *
  141.       IF( N.EQ.1 ) THEN
  142.          IF( ICOMPZ.EQ.2 )
  143.      $      Z( 1, 1 ) = CONE
  144.          RETURN
  145.       END IF
  146. *
  147. *     Determine the unit roundoff and over/underflow thresholds.
  148. *
  149.       EPS = DLAMCH( 'E' )
  150.       EPS2 = EPS**2
  151.       SAFMIN = DLAMCH( 'S' )
  152.       SAFMAX = ONE / SAFMIN
  153.       SSFMAX = SQRT( SAFMAX ) / THREE
  154.       SSFMIN = SQRT( SAFMIN ) / EPS2
  155. *
  156. *     Compute the eigenvalues and eigenvectors of the tridiagonal
  157. *     matrix.
  158. *
  159.       IF( ICOMPZ.EQ.2 )
  160.      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
  161. *
  162.       NMAXIT = N*MAXIT
  163.       JTOT = 0
  164. *
  165. *     Determine where the matrix splits and choose QL or QR iteration
  166. *     for each block, according to whether top or bottom diagonal
  167. *     element is smaller.
  168. *
  169.       L1 = 1
  170.       NM1 = N - 1
  171. *
  172.    10 CONTINUE
  173.       IF( L1.GT.N )
  174.      $   GO TO 160
  175.       IF( L1.GT.1 )
  176.      $   E( L1-1 ) = ZERO
  177.       IF( L1.LE.NM1 ) THEN
  178.          DO 20 M = L1, NM1
  179.             TST = ABS( E( M ) )
  180.             IF( TST.EQ.ZERO )
  181.      $         GO TO 30
  182.             IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
  183.      $          1 ) ) ) )*EPS ) THEN
  184.                E( M ) = ZERO
  185.                GO TO 30
  186.             END IF
  187.    20    CONTINUE
  188.       END IF
  189.       M = N
  190. *
  191.    30 CONTINUE
  192.       L = L1
  193.       LSV = L
  194.       LEND = M
  195.       LENDSV = LEND
  196.       L1 = M + 1
  197.       IF( LEND.EQ.L )
  198.      $   GO TO 10
  199. *
  200. *     Scale submatrix in rows and columns L to LEND
  201. *
  202.       ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) )
  203.       ISCALE = 0
  204.       IF( ANORM.EQ.ZERO )
  205.      $   GO TO 10
  206.       IF( ANORM.GT.SSFMAX ) THEN
  207.          ISCALE = 1
  208.          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
  209.      $                INFO )
  210.          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
  211.      $                INFO )
  212.       ELSE IF( ANORM.LT.SSFMIN ) THEN
  213.          ISCALE = 2
  214.          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
  215.      $                INFO )
  216.          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
  217.      $                INFO )
  218.       END IF
  219. *
  220. *     Choose between QL and QR iteration
  221. *
  222.       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
  223.          LEND = LSV
  224.          L = LENDSV
  225.       END IF
  226. *
  227.       IF( LEND.GT.L ) THEN
  228. *
  229. *        QL Iteration
  230. *
  231. *        Look for small subdiagonal element.
  232. *
  233.    40    CONTINUE
  234.          IF( L.NE.LEND ) THEN
  235.             LENDM1 = LEND - 1
  236.             DO 50 M = L, LENDM1
  237.                TST = ABS( E( M ) )**2
  238.                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
  239.      $             SAFMIN )GO TO 60
  240.    50       CONTINUE
  241.          END IF
  242. *
  243.          M = LEND
  244. *
  245.    60    CONTINUE
  246.          IF( M.LT.LEND )
  247.      $      E( M ) = ZERO
  248.          P = D( L )
  249.          IF( M.EQ.L )
  250.      $      GO TO 80
  251. *
  252. *        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
  253. *        to compute its eigensystem.
  254. *
  255.          IF( M.EQ.L+1 ) THEN
  256.             IF( ICOMPZ.GT.0 ) THEN
  257.                CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
  258.                WORK( L ) = C
  259.                WORK( N-1+L ) = S
  260.                CALL ZLASR( 'R', 'V', 'B', N, 2, WORK( L ),
  261.      $                     WORK( N-1+L ), Z( 1, L ), LDZ )
  262.             ELSE
  263.                CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 )
  264.             END IF
  265.             D( L ) = RT1
  266.             D( L+1 ) = RT2
  267.             E( L ) = ZERO
  268.             L = L + 2
  269.             IF( L.LE.LEND )
  270.      $         GO TO 40
  271.             GO TO 140
  272.          END IF
  273. *
  274.          IF( JTOT.EQ.NMAXIT )
  275.      $      GO TO 140
  276.          JTOT = JTOT + 1
  277. *
  278. *        Form shift.
  279. *
  280.          G = ( D( L+1 )-P ) / ( TWO*E( L ) )
  281.          R = DLAPY2( G, ONE )
  282.          G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
  283. *
  284.          S = ONE
  285.          C = ONE
  286.          P = ZERO
  287. *
  288. *        Inner loop
  289. *
  290.          MM1 = M - 1
  291.          DO 70 I = MM1, L, -1
  292.             F = S*E( I )
  293.             B = C*E( I )
  294.             CALL DLARTG( G, F, C, S, R )
  295.             IF( I.NE.M-1 )
  296.      $         E( I+1 ) = R
  297.             G = D( I+1 ) - P
  298.             R = ( D( I )-G )*S + TWO*C*B
  299.             P = S*R
  300.             D( I+1 ) = G + P
  301.             G = C*R - B
  302. *
  303. *           If eigenvectors are desired, then save rotations.
  304. *
  305.             IF( ICOMPZ.GT.0 ) THEN
  306.                WORK( I ) = C
  307.                WORK( N-1+I ) = -S
  308.             END IF
  309. *
  310.    70    CONTINUE
  311. *
  312. *        If eigenvectors are desired, then apply saved rotations.
  313. *
  314.          IF( ICOMPZ.GT.0 ) THEN
  315.             MM = M - L + 1
  316.             CALL ZLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ),
  317.      $                  Z( 1, L ), LDZ )
  318.          END IF
  319. *
  320.          D( L ) = D( L ) - P
  321.          E( L ) = G
  322.          GO TO 40
  323. *
  324. *        Eigenvalue found.
  325. *
  326.    80    CONTINUE
  327.          D( L ) = P
  328. *
  329.          L = L + 1
  330.          IF( L.LE.LEND )
  331.      $      GO TO 40
  332.          GO TO 140
  333. *
  334.       ELSE
  335. *
  336. *        QR Iteration
  337. *
  338. *        Look for small superdiagonal element.
  339. *
  340.    90    CONTINUE
  341.          IF( L.NE.LEND ) THEN
  342.             LENDP1 = LEND + 1
  343.             DO 100 M = L, LENDP1, -1
  344.                TST = ABS( E( M-1 ) )**2
  345.                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
  346.      $             SAFMIN )GO TO 110
  347.   100       CONTINUE
  348.          END IF
  349. *
  350.          M = LEND
  351. *
  352.   110    CONTINUE
  353.          IF( M.GT.LEND )
  354.      $      E( M-1 ) = ZERO
  355.          P = D( L )
  356.          IF( M.EQ.L )
  357.      $      GO TO 130
  358. *
  359. *        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
  360. *        to compute its eigensystem.
  361. *
  362.          IF( M.EQ.L-1 ) THEN
  363.             IF( ICOMPZ.GT.0 ) THEN
  364.                CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
  365.                WORK( M ) = C
  366.                WORK( N-1+M ) = S
  367.                CALL ZLASR( 'R', 'V', 'F', N, 2, WORK( M ),
  368.      $                     WORK( N-1+M ), Z( 1, L-1 ), LDZ )
  369.             ELSE
  370.                CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
  371.             END IF
  372.             D( L-1 ) = RT1
  373.             D( L ) = RT2
  374.             E( L-1 ) = ZERO
  375.             L = L - 2
  376.             IF( L.GE.LEND )
  377.      $         GO TO 90
  378.             GO TO 140
  379.          END IF
  380. *
  381.          IF( JTOT.EQ.NMAXIT )
  382.      $      GO TO 140
  383.          JTOT = JTOT + 1
  384. *
  385. *        Form shift.
  386. *
  387.          G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
  388.          R = DLAPY2( G, ONE )
  389.          G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
  390. *
  391.          S = ONE
  392.          C = ONE
  393.          P = ZERO
  394. *
  395. *        Inner loop
  396. *
  397.          LM1 = L - 1
  398.          DO 120 I = M, LM1
  399.             F = S*E( I )
  400.             B = C*E( I )
  401.             CALL DLARTG( G, F, C, S, R )
  402.             IF( I.NE.M )
  403.      $         E( I-1 ) = R
  404.             G = D( I ) - P
  405.             R = ( D( I+1 )-G )*S + TWO*C*B
  406.             P = S*R
  407.             D( I ) = G + P
  408.             G = C*R - B
  409. *
  410. *           If eigenvectors are desired, then save rotations.
  411. *
  412.             IF( ICOMPZ.GT.0 ) THEN
  413.                WORK( I ) = C
  414.                WORK( N-1+I ) = S
  415.             END IF
  416. *
  417.   120    CONTINUE
  418. *
  419. *        If eigenvectors are desired, then apply saved rotations.
  420. *
  421.          IF( ICOMPZ.GT.0 ) THEN
  422.             MM = L - M + 1
  423.             CALL ZLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ),
  424.      $                  Z( 1, M ), LDZ )
  425.          END IF
  426. *
  427.          D( L ) = D( L ) - P
  428.          E( LM1 ) = G
  429.          GO TO 90
  430. *
  431. *        Eigenvalue found.
  432. *
  433.   130    CONTINUE
  434.          D( L ) = P
  435. *
  436.          L = L - 1
  437.          IF( L.GE.LEND )
  438.      $      GO TO 90
  439.          GO TO 140
  440. *
  441.       END IF
  442. *
  443. *     Undo scaling if necessary
  444. *
  445.   140 CONTINUE
  446.       IF( ISCALE.EQ.1 ) THEN
  447.          CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
  448.      $                D( LSV ), N, INFO )
  449.          CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
  450.      $                N, INFO )
  451.       ELSE IF( ISCALE.EQ.2 ) THEN
  452.          CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
  453.      $                D( LSV ), N, INFO )
  454.          CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
  455.      $                N, INFO )
  456.       END IF
  457. *
  458. *     Check for no convergence to an eigenvalue after a total
  459. *     of N*MAXIT iterations.
  460. *
  461.       IF( JTOT.EQ.NMAXIT ) THEN
  462.          DO 150 I = 1, N - 1
  463.             IF( E( I ).NE.ZERO )
  464.      $         INFO = INFO + 1
  465.   150    CONTINUE
  466.          RETURN
  467.       END IF
  468.       GO TO 10
  469. *
  470. *     Order eigenvalues and eigenvectors.
  471. *
  472.   160 CONTINUE
  473.       IF( ICOMPZ.EQ.0 ) THEN
  474. *
  475. *        Use Quick Sort
  476. *
  477.          CALL DLASRT( 'I', N, D, INFO )
  478. *
  479.       ELSE
  480. *
  481. *        Use Selection Sort to minimize swaps of eigenvectors
  482. *
  483.          DO 180 II = 2, N
  484.             I = II - 1
  485.             K = I
  486.             P = D( I )
  487.             DO 170 J = II, N
  488.                IF( D( J ).LT.P ) THEN
  489.                   K = J
  490.                   P = D( J )
  491.                END IF
  492.   170       CONTINUE
  493.             IF( K.NE.I ) THEN
  494.                D( K ) = D( I )
  495.                D( I ) = P
  496.                CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
  497.             END IF
  498.   180    CONTINUE
  499.       END IF
  500.       RETURN
  501. *
  502. *     End of ZSTEQR
  503. *
  504.       END
  505.